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Abstract 

For the last almost three decades, since the famous Buchberger-Moller(BM) 
algorithm emerged, there has been wide interest in vanishing ideals of points 
and associated interpolation polynomials. Our paradigm is based on the 
theory of bivariate polynomial interpolation on cartesian point sets that gives 
us related degree reducing interpolation monomial and Newton bases directly. 
Since the bases are involved in the computation process as well as contained 
in the final output of BM algorithm, our paradigm obviously simplifies the 
computation and accelerates the BM process. The experiments show that 
the paradigm is best suited for the computation over finite prime fields that 
have many applications. 
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1. Introduction 

For an arbitrary field F, we let ¥ q a finite prime field of size q and U d := 
¥[x\, . . . , Xd] the (i-variate polynomial ring over F. Given a preassigned set of 
distinct affine points S C ¥ d , it is well-known that the set of all polynomials 
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in U d vanishing at S constitutes a radical zero- dimensional ideal, denoted by 
X(E), which is called the vanishing ideal of H. 

Recent years, there has been considerable interest in vanishingideals of 
points in many branches of mathematics such as algebraicgeometry[l|, multi- 
variate interpolation [i], Ijjj], coding theory (H, statistics p], and even compu- 
tational molecular biology 0, As is well known, the most significant mile- 
stone of the computation of vanishing ideals is the algorithm presented in [i[ 
by Hans Michael Moller and Bruno Buchberger known as Buchberger-Moller 
algorithm(BM algorithm for short). For any point set H C ¥ d and fixed term 
order -<, BM algorithm yields the reduced Grobner basis for X(S) w.r.t. -< 
and a -<-degree reducing interpolation Newton basis for <i-variate Lagrange 
interpolation on H. The algorithm also produces the Grobner escalier of X(E) 
w.r.t. -< as a byproduct. Afterwards, in 1993, BM algorithm was applied in 
To| in order to solve the renowned FGLM-problem. In the same year, 



merged BM and FGLM algorithms into four variations that can solve more 
general zero-dimensional ideals therefore related ideal interpolation problems 
[3j. The algorithms are referred as MMM algorithms. 

Although very important, BM algorithm (and MMM algorithms) has a 
very poor complexity that limits its applications. In this decade, many au- 
thors proposed new algorithms that can reduce the complexity but mostly 



suitable for special cases. [12| presented a modular version of BM algorithm 



that is best suited to the computation over Q. [13j, [14|, [15| presented algo- 
rithms for obtaining, with relatively little effort, the Grobner escalier of a 
vanishing ideal w.r.t. the (inverse) lexicographic order that can lead to an 
interpolation Newton basis or the reduced Grobner basis for the vanishing 
ideal after solving a linear system. 

For a fixed point set S in ¥ d and a term order -<, it is well known that 
there are two factors that determine the Grobner escalier of X(E) w.r.t. -< 
thereby the reduced Grobner basis for X(S) and related degree reducing in- 
terpolation Newton bases (up to coefficients). One is apparently the cardinal 
of H. It is the unique determinate factor in univariate cases. Another one is 
the geometry (the distribution of the points) of S that is dominating in multi- 
variate cases but not taken into consideration by BM and MMM algorithms. 
Recent years, [l6|, 17, 18] studied multivariate Lagrange interpolation on a 
special kind of point sets, cartesian point sets (aka lower point sets), and 
constructed the associated Grobner escalier and degree reducing interpola- 
tion Newton bases theoretically. We know from (9), [ll[ that, for a cartesian 
subset of 5 (it always exists!), certain associated degree reducing interpola- 
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tion Newton basis forms part of the output of BM algorithm w.r.t. some 
reordering of S. Therefore, finding a large enough cartesian subset of S with 
little enough effort will reduce the complexity of BM algorithm. 

Following this idea, the paper proposes a preprocessing paradigm for BM 
algorithm with the organization as follows. The next section is devoted as a 
preparation for the paper. And then, the main results of us are presented in 
two sections. Section [3] will pursue the paradigm for two special term orders 
while Section H] will set forth our solution for other more general cases. In the 
last section, Section [5j some implementation issues and experimental results 
will be illustrated. 



2. Preliminary 

In this section, we will introduce some notation and recall some basic 
facts for the reader's convenience. For more details, we refer the reader to 
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We let No denote the monoid of nonnegative integers. A polynomial f £ 
II 2 is of the form 

f=J2 f<* X "> #{« e K : + f a £ F} < oo, 

where monomial X a = x ai y a ' 2 with ct = (a^a^). The set of bivariate 
monomials in II 2 is denoted by T 2 . 

Fix a term order -< on IT 2 that may be lexicographical order -<i cx , in- 
verse lexicographical order -<i n i GX , or total degree inverse lexicographical order 
-<tdinicx etc. For all / £ II 2 , with / ^ 0, we may write 

/ = f^X* + f l2 X^ + ■■■ + f^, 

where ^ / Ti £ F, 7* £ Ng, i = 1, . . . , r, and X^ 1 y > y X^. We 

shall call LT(/) := / 7l X 71 the leading term and LM(/) := X 11 the leading 
monomial of /. Furthermore, for a non-empty subset F C II 2 , put 

LT(F) := {LT(/) : / £ F}. 



As in 21[, we define the -<— degree of a polynomial / £ II 2 to be the leading 
bidegree w.r.t. -< 

«y(/):=7, X~< = LM(f), 
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with 5(0) undefined. Further, for any finite dimensional subset F C II 2 , 
define 

5(F) :=max5(f). 

Finally, for any /, g G II 2 , if 5(f) -< 5(g) then we say that / is of lower degree 
than g and use the abbreviation 

f*g:= 5(f) -< 5(g). 

In addition, / ^ g is interpreted as the degree of / is lower than or equal to 
that of g. 

Let A be a finite subset of Ng. A is called a lower set if, for any ex. = 
(ai, 0:2) £ A, we always have 

R(a) := a 2 ) G Nq : < «■ < a i; z = 1, 2} C A. 

Especially, G A. Moreover, we set rrij = maxrh^eA h, < j < v, with v = 
max( 0i fe)(=^ A;. Clearly, A can be determined uniquely by the ordered (v + 1)- 
tuple (m , mi, . . . , m„) hence represented as L a ,(m , mi, ... , m v ). Swapping 
the roles of x and y, we can also represent A as L y (n , m, n mo ) with 
rij = max(^) e< 4 /c, < i < m . It should be noticed that v = no. 

Given a set S = . . . , ^ < - A1 - ) } C F 2 of fx distinct points. For prescribed 
values fi G F, i = 1, . . . , fx, find all polynomials p G II 2 satisfying 

p(£ (i) ) = /i, i = 1,-.-,^- (1) 

We call it the problem of bivariate Lagrange interpolation. Note that in most 
cases, especially from a numerical point of view, we are not interested in all 
such p's but a "degree reducing" one, as in the univariate cases. 

Definition 1. 0] Fix term order -<. We call a subspace V C II 2 a degree 
reducing interpolation space w.r.t. -< for the bivariate Lagrange interpolation 
©if 

DR1. V is an interpolation space, i.e., for any /j G F, i — 1, . . . , fx, there is a 
unique p G V such that p satisfies (JT]). In other words, the interpolation 
problem is regular w.r.t. P. 

DR2. P is -<— reducing, i.e., when denotes the Lagrange projector with 
range V, then the interpolation polynomial 

L v q <q, Wq G n 2 . 
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For interpolation problem flTJ, a given interpolation space V C II 2 will 



give rise to an interpolation scheme that is referred as (5, V), cf. [20|. Since 
(JTJ) is regular w.r.t. P, we can also say that (S, P) is regular. Moreover, if T 7 
is degree reducing w.r.t. -<, a basis {pi, . . . ,p^} for V will be called a degree 
reducing interpolation basis w.r.t. -< for ([T]). Assume that p\ -< P2 -< ■ ■ ■ -< p M . 
If 

P;(£ W ) = %, 1 < * < i < A*, 

for some suitable reordering of S, then we call {pi, . . . , p^} a degree reducing 
interpolation Newton basis(DKINB) w.r.t. -< for ([T|). 

Let be the reduced Grobner basis for the vanishing ideal X(E) w.r.t. -<. 
The set 

N^(X(£)) := {x" G T 2 : LT(</) f G G^} 

is called the Grobner escalier of X(H) w.r.t. -<. From (2!, Hi], the interpolation 
space spanned by N^(X(S)), denoted by is canonical since it is the 

unique degree reducing interpolation space spanned by monomials w.r.t. -< 
for (JTJ) . Hence, we call N^(X(S)) the degree reducing interpolation monomial 
basis(DRlMB) w.r.t. -< for with #N^(X(S)) = p. Let 

N X (S) := {ex. : x a G N^(X(S))} C N 2 . 

We can deduce easily that (S) is a lower set and obviously has a one-to-one 
correspondence with N^(X(S)). Therefore, interpolation scheme (S,"P^(S)) 
can be equivalently represented as (H, N_<(5)). 

According to [17|, we can construct two particular lower sets from S, 
denoted by S X (E), S y (E), which reflect the geometry of H in certain sense. 

Specifically, we cover the points in S by lines 1^,1*, ... ,1* parallel to the 
x-axis and assume that, without loss of generality, there are rrij+1 points, say 
Uqj, ufj, . . . , Um jt j, on I* with m > mi > ■ ■ ■ > m v > hence the ordinates 
of ufj and u x V pi ^ i', same. Now, we set 

S X (E) :={(t,j):0<t<m J , < j < u}, 

which apparently equals to L x (m , mi, . . . , m u ). We can also cover the points 
by lines 1%, If, . . . , l v x parallel to the y-ecxis and denote the points on line l\ by 
uf , u^, ... , ixf n . with n >«!>■■■> ra A > hence the abscissae of u\a and 
w^,, j 7^ j', same. Similarly, we put 

S y (E) := {(i,j) : < i < A, < j < m} = L„(n ,ni, . . . ,n A ). 
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In addition, we can also define the sets of abscissae and ordinates 

H j (E):={x:(x,y)el*nE}, < j < u, 
K(H):={j/:(5,j/)GZj'nH}, < z < A. 



Definition 2. |17| We say that a set E of distinct points in F 2 is cartesian 
if there exists a lower set A such that S can be written as 

S = {(xi,yj) : (i,j) G A}, 

where the distinct numbers, and similarly the y/s. We also say that 

E is ^.-cartesian. 

To the best of our knowledge, there are two criteria for determining 
whether a 2-dimensional point set is cartesian. 

Theorem 1. [TtJ A set of distinct points S C F 2 is cartesian if and only if 

S X (~) = Sy(Z). 



Theorem 2. [18j -A set o/ distinct points S C F 2 cartesian if and only if 
H (E) D ffi(S) D • • • D i7,(H), 7 (S) D ^i(H) D • • • D 7 A (S). 



About the bivariate Lagrange interpolation on a cartesian set, [17| proved 
the succeeding theorem. 



Theorem 3. 17| Given a cartesian set S C F 2 , there exists a unique lower 
set A G Nq such that S zs A- cartesian and the Lagrange interpolation scheme 
(S, ^4) zs regular. 

Finally, we will redescribe the classical BM algorithm with the notation 
established above. 

Algorithm 1. (BM Algorithm) 

Input: A set of distinct points S = : % = 1, . . . , fi} C ¥ d and a fixed 
term order -<. 

Output: The 3-tuple (G,N,Q), where G is the reduced Grobner basis 
for 1(E) w.r.t. -<, N is the Grobner escalier of 1(E) (the DRIMB for $T]) 
also) w.r.t. -<, and Q is a DRINB w.r.t. -< for 

BM1. Start with lists G = [ ],N = [ ],Q = [ ],L = [1], and a matrix 
I? = (b^) over F with \i columns and zero rows initially. 
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BM2. If L = [ ], return (G,N,Q) and stop. Otherwise, choose the 
monomial t = min^L, and delete t from L. 

BM3. Compute the evaluation vector . . . , and reduce it 

against the rows of B to obtain 

(v 1 ,...,v lt ) = • • • , t(^)) - <b a , 6i M ), a* G F. 

i 

BM4.. If (v\, . . . , u M ) = (0, . . . , 0), then append the polynomial t—^2 i 
to the list G, where is the zth element of Q. Remove from L all the multiples 
of t. Continue with BM2. 

BM5. Otherwise . . . , v^) ^ (0, . . . , 0), add (v±, . . . , v^) as a new row 
to B and t — Y2i a iQi as a new element to Q. Append the monomial t to N, 
and add to L those elements of {xit, . . . , x^t} that are neither multiples of 
an element of L nor of LT(G). Continue with BM2. 

3. Special cases 

In this section, we will focus on -<i cx and -<i n i cx that may be the most 
talked about term orders. For these special cases, our preprocessing paradigm 
will first provide exact N, Q of the 3-tuple output (G, N, Q) to BM algorithm 
directly and effortlessly. And then, G can be obtained by BM algorithm 
easily. Note that we will continue with all the notation that we established 
for S X (E) and S y (E) in the previous section. 

Proposition 4. Let H be a set of n distinct points u x mn = (x mn ,y mn ) G 
F 2 , (m,n) G S X (E). The points give rise to polynomials 

i-i i-l 
t=0 s=0 

where ip^ = 1/ YllZoiVoj ~ Hot) YY s Zo( x ij ~ x sj) e F > an d the empty products 
are taken as 1. Then we have 

tfjiUmn) = %i),(m,n), (h j) ^inlcx (m,n). 

Proof. Fix G S X (E). Recalling the definition of u^, we have yoj = yij. 
If = (m, n), by y o 7^ 2/oi 7^ ■ ■ ' 7^ Uoj and ^oj 7^ x ij 7^ • • • 7^ x ij, we nave 

i-l i-l i-l 

s=0 t=0 s=0 



t=0 
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which implies 4>ij( u ij) — 1- 

Otherwise, if ^inicx ( m ,n), we have j > n, or j = n, i > m. When 
j > n, we have 



i-l 

^■(«m„) = VijiVmn ~ Voo) ' " ' (j/mn ~ 2A)n) " ' ' (l/mn - VOJ-l) \\{x mn - X sj ) 

i-l 

= VijiVOn ~ VOO) ■ ■ ■ {VOn ~ VOn) • ■ ■ (VOn ~ UOj-l) JJOmn - X sj ) 

= o, 

and when j = n,i > m, 

i-i 

<i>ij«n) = ^jY[(ymn 
t=0 

n-1 

= <Pij WiVmn 
t=0 

= 0, 

which leads to 

4>ij{Kin) = °> (hj) ^inlox (m,7l). 

□ 

Similarly, we can prove the following proposition: 

Proposition 5. Let H be a set of \i distinct points u y mn = (x mn ,y mn ) G 
F 2 , (m,n) G S y (3). We define the polynomials 

i-i i-i 

s=0 t=0 

w/iere <^ = 1/ ^=0(^0 - ^so) l~It=o(^i ~ Hit) e F - r/ie em Pty products are 
taken as 1. Then, 

$j«m) = %),(m,n), ^lex (m,7l). 



yot)^mn ^Oj) ' ' ' y^mn •^mj) ' ' ' (^mn l,i) 
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In 2004, [17]] proved that the Lagrange interpolation schemes (E,S X (E)) 
and (S, S y (E,)) are both regular. Here we reprove the regularities in another 
way for the purpose of presenting the degree reducing interpolation bases 
theoretically . 

Theorem 6. Resume the notation in Proposition 0] and O Then the La- 
grange interpolation schemes (H, S X (E)) and (H, 5 y (H)) are regular. Further- 
more, 

(i) the setN x := {x i y j : (i,j) G S X (E)} is the DRIMB as well as Q x := {<%• : 

G S X (E)} is a DRINB w.r.t. -<i cx for the interpolation problem (Tj[|). 

(ii) the set N y := {x l y^ : G S y (H)} t/ie DRIMB as well as Q y : = 
{4> v i3 : G S„(S)} is a DRINB w.r.t. ^ inlcx for ©. 

Proof. We only give the proof for S X (E). The statements about S y (E) can 
be proved likewise. 

First, we will show the regularity of the interpolation scheme (S, S X (E)). 
Let V x := Spa%iV x C II 2 with dimP^, = #S = fi. Obviously, N x is the 
monomial basis for it. By (j3J), we can check easily that 

SpanpQz C V x . 

Construct a square matrix B^ XfJi whose (h,k) entry is </>h( u fc) where 4>hi u t 
are hth and kth elements of Q x and H = {w^n : ( m , n ) G S^S)} w.r.t. the 
increasing -<i n i cx on (z,j) and (m, n) respectively. From Proposition HI B^ XfJi 
is upper unitriangular which implies that Span F Qa. = V x and Q x forms a 
Newton basis for V x . It follows that V x is an interpolation space for Lagrange 
interpolation (JT]) therefore the scheme (E,V X ) is regular. Since (E,S X (E)) = 
(E,V X ), according to Section [21 (H, S X (E)) is regular. 

Next, we shall verify that the statements in (i), which is equivalent to the 
statement that V x is a degree reducing interpolation space w.r.t. -<i cx for (Tj[|) 
that coincides with P^ lox (S). Since the arguments above have proved that 
V x satisfies the DR1 condition in Definition [H what is left for us is to check 
the DR2 condition. From j^H, we only need to check it for monomials. 

Take a monomial x io yi° G T 2 . We shall prove that 

L-p x x io y jo ^icx zV°- (5) 

Since V x satisfies DR1, L-p x x l °y^° is the unique polynomial in V x that matches 
x^yio on 2. Therefore, when x l °y^° G N x , we have L-p x x l °yi° = x lo y 30 , 
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namely (jSJ) is true for this case. Assume that 

S x {^) — L x (mo, • • • , rn no ) = L y (no, . . . , n mo ). 

It is easy to see that 5{V X ) = (m ,n mo ). If x m °y nm o -< lex x io y jo then 
5(L Tx x io y jo ) ^i cx 6{P X ) = {m ,n mo ) -< lex (i ,jo) = 5(x io y jo ) that leads to 
(JSJ) for the case. 

Thus, what remains for us is to check (JSJ) for x lo yi° ^ with (io, Jo) ~<icx 
(mo, n mo ) that implies < io < m ,j > n io . For this, we only need to verify 
that 

L Vx x io y j ° G Span F {xy : G F, }, (6) 

where F l0 = G S X (E) : ^ lcx (z ,Jo)} C S X (E). If xV° G 1(E), 

then L-p x x to y^ = -<i ex x l °y^°. The statement fl6]) becomes trivial in this 
case. Otherwise, if we can find a polynomial p G II 2 such that 

p = x io y io_ Oijx'y* el(E), (7) 

where a^- G F are not all zero, then ([6]) follows. 

According to Section [2j our point set E = {u^ = (xij,yij) : G 
S X (E)}. Let E' = {u x mn G H : (m,n) G F io } c S. Now, we claim that there 
exists a unique polynomial p of the form (J7|) such that p G X(E'), which is 
equivalent to the statement that the linear system 

Ea- V v j = r io v jo v x G CSl 
^13 mni)mn ^mntlmni ""ran *~ '— ' ' V"/ 

(«J)6^0 

has a unique solution. 

Note that Span^x 4 ?/- 7 : G F io } = Span F {0f- : G i 7 ^}. We can 

conclude that the rank of the coefficient matrix of ([H]) is equal of that of the 
matrix BL F . X _^ F . , which is a submatrix of B whose (h,k) entry is 4>%{ u k) 
where <f>%,u% are hth and kth elements of {0?- : G i<i } and E' = {u^} 
w.r.t. the increasing -<i n i cx on and (m, n) respectively. By (J3J), we see 
easily that B' is upper unitriangular which implies that the coefficient matrix 
of dSJ) is of full rank. Accordingly, there is a unique polynomial p G X(E') 
that has the form (J7J). 
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Now we shall verify that p(u^) = 0, G H \ H'. By the definition of S', 
we know that i > i$ here. Let 

q {x) := p(x, yij ) = b ^ s e nl > Ke¥. 

Since y 0j = yij = ■ ■ ■ = Vi ,j = Vij and u^, u x j, . . . , u x ia j G E', it follows that 

<?Osj) = P(x s j, Vij) = p(x sj ,y sj ) = p(u x sj ) = 0, s = 0, . . . , z' , 

namely has io + 1 zero points that clearly implies = 0. Since 
P( u ij) = l( x ij) = 0' we have p G 1(5). By (EJ), ((5]) is true in this case. As a 
result, for any / G II 2 , we have 

that is to say V x satisfies DR2. 

Consequently, by Definition [1], V x is a degree reducing interpolation space 
w.r.t. -<i cx for Lagrange interpolation ([T|). Hence N x is the DRIMB and Q x 
is a Newton basis w.r.t. -<i cx for (FTj). □ 



Note that "P_< lex (H) is the unique degree reducing interpolation space 
spanned by monomials w.r.t. -<i CX! thus we have V x = P^ lcx (H). Therefore, 
N x = A r ^ lox (X(S)) holds, which means that N x is also the Grobner escalier of 
1(E) w.r.t. -<i ex . 

Corollary 7. IfEc F 2 is an A-cartesian set, then A = S X (E) = S y (E). 

Proof. Since 5 is cartesian, by Theorem [1] and El we have S X (E) = S y (E) 
hence (E,S X (E)) = (S, S y (E)) are both regular. But from Theorem [3j only 
A can make (E, A) regular, therefore A = S X (E) = S y (E). □ 

From Algorithm [T] we know that G, N, Q are essential elements of BM 
algorithm and compose its output. For -<i ex and -<i n i cx cases, Theorem [6] 
presents us N and Q theoretically hence we can obtain them with little 
effort. According to 11], the leading terms of G are contained in the border 
set of N. Therefore, we can get G faster than compute G directly with BM 
algorithm. Now is our algorithm. 
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Algorithm 2. (SPBM) 

Input: A set of distinct affine points 3 C F 2 and fixed -<i cx or -<i n i ex . 

Output: The 3-tuple (G, N, Q), where G is the reduced Grobner basis of 
X(3), N is the Grobner escalier N(Z(S)), and Q is a DRINB for the Lagrange 
interpolation on 3. 

SPBM1. Construct lower set S X (E) or S y (E) according to Section [2J 

SPBM2. Compute the sets iV and Q by Theorem 

SPBM3. Construct the border set L := {x-t : t E N} \J{yt : t E N}\N 
and the matrix B that is same to the in the proof of Theorem [61 

SPBM4. Goto BM2 of BM algorithm for the reduced Grobner basis G. 

Example 1. Let 

3 = {(0,1), (0,3), (1,0), (1,2), (1,3), (1,4), (2,1), (2,2), (3, 1)} C Q 2 . 

First, we choose lines x = 1, x — 0, x = 2, x = 3 as respectively 
(Shown in (a) of Figured]), therefore we have 

S y = {(0, 0), (0, 1), (0, 2), (0, 3), (1, 0), (1, 1), (2, 0), (2, 1), (3, 0)}, 

which is illutrated in (b) of Figure [TJ 



o i\ 



'02 



I 
I 
I 



I I 



♦- 

■ i 



(a) S (b) S y 

Figure 1: The point set and related S y of Example [TJ 
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Thus, by Theorem 0, we have 



N={l,y,y ,y ,x,xy,x , x y, x }; 

„ n 1 1 2 2 1 3 5 2 3 n 1 1 1 1 

Q — il,—v,—y v,—v V H — V, — £ + 1, — xy H — yH — a; . 

i ' 2 y '3 y 3 8 8 4 2 2 2 2' 

1 2 1 1 2 1 1 2,1 1 s 1 2,1 x 
-x x,-x y xy x H — x, -x x H — x . 

2 2 ' 2 y 2 y 2 2 ' 6 2 3 J 



Next, from SPBM3, the border set L = {y 4 , xy 2 , xy 3 , x 2 y 2 , x 3 y , x 4 } and the 
matrix 

/ 1 1 1 ... \ 
1 3/2 ■•• 
1 ■•• 



B 



\ : : : -.J 

Finally, turn to BM2 with these N, Q, L, B, we can get the reduced Grobner 
basis 



G = {x 4 — 6x 3 + llx 2 — 6x, x 3 y — 3x 2 y + 2xy — x 3 + 3x 2 — 2x, 
xy 2 — y 2 + ~~ 7^ X V + 4y — -x 2 + — 3, 
y 4 - 9y 3 + 26y 2 - \ 2 y + —xy - 27y - 3x 3 + — x° 



yX + 9}. 



for 1(H) w.r.t. -< in i cx . 

Example 2. Given a bivariate point set 

2 = {(0, 0), (0, 2), (0, 3), (1, 1), (5 0), (5 1), (5 2), (4, 0), (4, 2)} C Q 2 . 

We choose lines y = 0, y = 2, y = 1, y = 3 as l^jfj^^ respectively (Illus- 
trated in (a) of Figure |2]), which follows that 



S x = {(0, 0), (1, 0), (2, 0), (0, 1), (1, 1), (2, 1), (0, 2), (1, 2), (0, 3)}. 
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-f 



im; 



(a) S (b) ^ 

Figure 2: Illustrations for Example [21 



Thus, with SPBM algorithm, we have 



N 
Q 



n 2 2 2 2 3l 

{l,x,x ,y,xy,x y,y ,xy ,y }, 

1 4 2 16 1 1 2 2 8 2 

x H x,-y,-xy, x y H xy,—y +2y, 

4 15 15 ' 2 y '8 y 15 y 15 y y y 



2 q 2 o 4 4 1 q 

-xy + -y + -xy - -y, -y 

6 6 6 6 



-,y}, 



G = {y — Qy +lly — Qy, xy — 3xy + 2xy, x y —2xy — -xy 

5 o ^^2 1^ S o2 ^ 

+ <xy - -y + — y - — y, x - —x - 3xy + Qxy 



lOx 



15 o 75 2 45 , 

— y H y yj. 

4 4 2 J 



4. General cases 

Next, we will discuss how to accelerate BM algorithm with respect to 
term orders other than -<i cx or -<i n i cx - In [17J, the author proposed that if 
the set of points S is cartesian, then we can obtain the interpolation basis 
without any difficulty, see Theorem [3j But in general H may not be cartesian. 
However, we have the following proposition. 

Proposition 8. There must exist at least one cartesian subset for any non- 
empty set of points in F 2 . 
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Proof. Let S be a non-empty set of points. Hence, there exists at least one 
point £ G S. But £ itself can construct a cartesian subset {£} C H. □ 

Definition 3. Let 5 be a set of points in F 2 and H' be a cartesian subset of 
H. We say that H' is a maximal cartesian subset of S if any cartesian proper 
subset H" of H containing S' is such that H" = S'. In addition, a maximal 
row subset of S is a non-empty subset that equals the intersection of S and 
a horizontal line. 

From Proposition [8] we know that, for a set of given points, we can surely 
find a maximal cartesian subset of it. Is it unique? Unfortunately, the answer 
is often false. 

Example 3. Recall Example [2} let 

{(0, 0), (0, 2), (| 0), (| 1), (| 2), (4, 0), (4, 2)}, 

{(0,0),(0,2),(0,3),(|o),(|2),(4,0),(4,2)}, 

{(l,l),(|0),(|l),(|2)}. 

We can check easily that S'^S^Eg are all maximal cartesian subsets of S 
(Illustrated in Figure EJ. 

Lemma 9. Let S be a set of distinct points in F 2 and -< a fixed term order. 
IfE' is an A' -cartesian subset o/S, then 

i' = N,(H')cN^), 

or equivalently, 

{*y : G A'} = N^(Z(H')) C N^(Z(~)). 

Proof. From Section^ the Grobner escalier N^(X(H')) is the DRIMB w.r.t. 
-< for the bivariate Lagrange interpolation on H' hence the interpolation 
scheme (H', N^(S')) is regular. Since ^4' C Nq is lower and H' is v4'-cartesian, 
according to Theorem [31 A' is the unique lower set making the bivariate 
Lagrange interpolation on S' regular. This gives 

A' = N^(E'). 
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Figure 3: Maximal cartesian subsets of S, where • denotes the points in 
while o denotes the points in 
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Since S' C E, from [19j, we know that the vanishing ideals satisfy X(H') D 
X(S). Denote by G', G the reduced Grobner bases for X(S') and X(2) w.r.t. -< 
respectively. We will prove N^(X(E')) C N_<(X(S)) by contradiction. For any 
x % yi G N^(X(H')), we suppose there were some g e G such that UT{g)\x l y^ . 



By [19 



(LT(G')> = (LT(X(S')) D LT(X(S)) D LT(G). 

Therefore, LT(#) e LT(G) C (LT(G')) implies that there exists some g' e G' 
such that LT(^')|LT(gf). Since \JY{g)\x l y\ we have LT (g^^yi that contra- 
dicts our assumption on x l y\ which proves that N_<(X(S')) C N X (X(H)) due 
to the definition of N X (X(S)). Finally, N X (H') ^ N^(X(S')) and N^(S) = 
N^(X(S)) complete the proof. □ 

Remark 1. For any ^4-cartesian set E, by CorollaryO we have A = S X (E) = 
S y (E) that obviously leads to A= S X (E) = S y (E) = N_<(H), according to the 
Lemma above, where term order -< is arbitrary. 

Now comes an algorithm for constructing a maximal cartesian subset of 
a given point set in F 2 . 

Algorithm 3. (Maximal Cartesian Subset Construction Algorithm) 

Input: A set of distinct points E = {£W : z = 1, . . . , //} c F 2 . 
Output: A maximal cartesian subset E' of E. 
MCS1. Start with an empty list E' = [ ]. 

MCS2. If E = [ ], return the set S' and stop. Otherwise, compute lower 
sets S X (E) and S y (E). 

MCS3. If S X (E) = S y (E), then replace E' by S' U S, return the set 5' 
and stop. 

MCS4. Otherwise, we first choose a maximal row subset of E with 
maximal cardinal number, denoted by A. Next, delete from E the points 
either in A or have different abscissae from the points in A. Finally, replace 
E' by E' U A and continue with MCS2. 

The following theorem ensure that this algorithm will terminate in finite 
steps with a maximal cartesian subset as its output. 

Theorem 10. The algorithm described above will stop in a finite number 
of loops. Furthermore, the set E' returned by the algorithm is a maximal 
cartesian subset. 



17 



Proof. As input data of Algorithm [21 point set E is finite. Observing that 
#E decreases actually in every loop, the algorithm will terminate in a finite 
number, say M, of loops for sure. We assume that M > 1 since M — 1 is 
trivial. 

Ej n and S^ ut signify the input and output E' of MCS4 in some loop 
respectively. Next, we will prove by induction on 1 < r < M — 1 that in the 
rth loop E' out is a cartesian set. The case r = 1 is obvious since E' in = [ ] and 
E' ont is clearly cartesian as a maximal row subset of E. Assume the statement 
is true for r = £ < M — 1. When r = / + 1, by the induction hypothesis, S( n 
is cartesian. Therefore, by Corollary [3, we assume that 

S; n = {(xi,yj) : e S x (E' in )}, 

where S x (E[ n ) = L x (m , . . . ,m no ) = L y (n , . . . , n mo ). Observing the con- 
struction process of E' in the algorithm, we see easily that n = n\ = 
■ ■ ■ = n m . Let the maximal row subset of H we choose at this moment 
be A = {(x^°\y), (x^\y), . . . , (x^ k \y)}. Due to the nature of A, we have 
k < m no and y ^ y j} j = 0, . . . , n . 

We claim that the set E[ n U A is cartesian. In fact, we will focus on the 
horizontal parallel lines lj : y = yj, j = 0, . . . , n , and ^ 0+1 : y — y. Resume 
the notation in ([2]). ifj(S' n UA) = Hj(E' in ) = {xi : < i < rrij},j = 0, . . . , n , 
and if no+ i(S( n U A) = {afW : < 2 < &;}. Since S[ n is S' z (S[ n )-cartesian, by 
Theorem H the relation # (H- n U A) D #i(£( n U A) D • ■ • D H no (E' in U A) 
holds. From the description of MCS4, we can deduce that H no (E' in U A) D 
if no+ i(5| n U A), which leads to 

i/ (~; n UA)D Hi(E[ n U A) D • • O ^„ 0+1 (S; n U A). (9) 

Note that for any < i < k, there exists hi E {0, 1, . . . , m no ] such 

that = Xhi- Therefore, we could find a permutation a of {0, 1, . . . ,m } 
satisfying a(i) — hi, % — 0, . . . , k, and a(i) — i,i — m no + 1, . . . , m . Choose 
lines l\ : x = x a {i), i = 0, . . . , m , that give rise to Vi(S. n ) = {yj : < j < 
n a (i)},i = 0,...,m . Since n = ni = ■ ■ ■ = n m „ o , the relation Vo(H-J = 
v i(-m) = ■■■ = V mno (E[ n ) D V mno+l (Zr) D ■ ■ ■ D V mo (~[ n ) holds. Observing 
that Vi{E[ n U A) = Vi{BJ U {y}, i = 0,...,k, and Vi(H- n U A) = ^(SJJ, % = 
k + 1, . . . , m , it is easy to get 

V (E[ n U A) = ■ ■ ■ = V k (E[ n UA)D V k+1 (E{ n U A) D ■ ■ ■ D V rno (E[ n U A). 
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Thus together with (jUJ), E' out = 5- n U A is cartesian due to Theorem [21 hence 
our statement is true. 

For the Mth loop, if S = [ ], then S' here equals to the S^ ut of the MCS4 
of the (M— l)th loop that is cartesian due to the statement above. Otherwise, 
since the algorithm stops in MCS3 of this loop, H is a non-empty cartesian 
set. Similar to the arguments above, we can prove that S' = E' out U H is also 
cartesian. 

Finally, we should verify that the output S' of the algorithm is maximal. 
Otherwise, there must exist a maximal ^(S'^-cartesian subset H" of S sat- 
isfying H" g S'. Take a point £ = (x i(V y jo ) with (z' ,Jo) = min^ tolex {(i, j) G 
S X (E") : (xj, G S" \ H'}. Suppose there exists a point in S' sharing the or- 
dinate with £o- If h is chosen as a point in the maximal row subset in MCS4 
of some loop, by the definition of £o ; we know that £o is surely contained 
in the set H of that step, which contradicts the definition of the maximal 
row subset. Otherwise, it must appear in the cartesian set S in MCS3 in 
the final loop. Then, by the definition of £o> it should be contained in H 
hence the output set H', which introduces a contradiction. If there does not 
exist a point in 5' sharing the ordinate with £ , since S" is also cartesian, 
by Theorem [21 it is easily to see that £o must remain in S in every loop, 
which contradicts the termination condition. As a result, the output of the 
Algorithm [3] is a maximal cartesian subset. □ 

Let us continue with the setup and notation in Algorithm [31 and assume 
that the final output of it is S' who is S^S'J-cartesian . We now discuss how 
to preprocess the BM algorithm with the help of S'. 

Define an order -<= on the set S. Let G S. We say that £W -<= 

£( 2 ) if one of the following conditions holds: 

(1) G S', and e (2) G H\S'. 

(2) = (aci!,^),^ = (x i2 ,y j2 ) G S' and (kji) -< in i cx (12,^2) with 
(ifc.Jfc) eS x (E%k = l,2. 

It should be noticed that the order is not total. For the points in 5\S', any 
order of them can be interpreted as increasing. Hereafter, we will suppose 
that the points in S = . . . , C have been ordered increasingly w.r.t. 
-<=, namely ^ -<■= £^\0 < i < j < By the definition of -<=, we have 
H' = {£ (1) ,...,£ (#s,) }- 
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According to Lemma [HI N' = \x l yi : G S X (E')} C N, with N as a 
member of the 3-tuple output of BM algorithm. Thus the other monomials 
of N are obviously contained in T 2 \N'. Notice that the generators of T 2 \N' 
are located in the border of N', denoted by L, we can continue to spot the 
elements in L by BM algorithm to complete N. 

Next, we will pay attention to the computation of the Newton basis. Since 
H' is cartesian, recalling Proposition HJ we can construct the polynomials 0?- 
w.r.t. S x (Ef). Order G S x (Ef), increasingly w.r.t. under -<i n icx, 

and denote them as qi, #2, ■ ■ ■ , Q#s>- Set the matrix 



/ 



B 



? 2 (e (i) : 



?i(e (2) ) 
g 2 (e (2) ) 



9i(£ (# ^ 



\ 



(10) 



By Proposition HI B is obviously upper unitriangular which implies that the 
polynomials gi> 92, • • • , g#~y constitute a Newton basis for P^(H') = Span F A"'. 

All in all, with the notation above, we get our preprocessing procedure 
for BM algorithm. 

Algorithm 4. (GPBM) 

Input: A set of distinct points H C F 2 and a term order -<. 
Output: The 3-tuple (G,N,Q). 

GPBM1: Get a maximal cartesian subset H' of H by Algorithm [3J 
GPBM2: Compute the lower set S X (E') w.r.t. S', the set N := {x^- 7 ' : 
G ^(H')}, and the set Q := {91,92, • • • ,9#s'} where the g^'s are as in 

<m- 

GPBM3: Construct L := {x ■ t : t G A} |J{y ■ t : t £ N} \ N and the 
matrix i? that is same to ffTU]) . 

GPBM4: Goto BM2 of the BM algorithm to complete the computation 
and get the whole output. 



5. Implementation and Timings 

From the above section, we can see easily that our preprocessing paradigm 
is more suitable to the cases where the constructed maximal cartesian subset 
S' forms a relatively large proposition in S. Especially, when the field F is 
finite, our preprocessing will play a more important role in consideration of 
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the nature of finite fields. In this section, we will present some experimental 
results to compare the effectiveness of our paradigm with the classical BM. 
First see an example with point set of small size. 

Example 4. We choose the field F7, and let 

E = {(0,0), (0,1), (0,4), (0,5), (1,0), 
(1,1), (1,4), (1,6), (2,1), (2,2), 
(2, 6), (3, 2), (4, 2), (4, 5), (4, 6), 
(5,1), (5, 5), (5, 6), (6,0), (6, 2)}. 

By Algorithm [3J we can construct the maximal cartesian subset 

5' = {(0, 1), (1, 1), (2, 1), (5, 1), (1,6), (2,6), (5,6), (1,0), (1,4)} 



hence get 



N ={1, x, x 2 , x 3 , y, xy, x 2 y, y 2 }, 

Q ={1, x, Ax 2 + 3x, 2x 3 + x 2 + Ax, 3y + 4, 3xy + Ax + Ay + 3, 

2x 2 y + bx 2 + xy + 6x + Ay + 3, Qy 2 + 1, 2y 3 + by}, 

L ={y A , xy 2 , xy 3 , x 2 y 2 , x 3 y, x A }, 

( 1 1 1 ••■ \ 
1 2 ■■• 



B 



1 ••• 

v j 

Put these N, Q, L, B into BM algorithm, we can get the final output 

N ={l,x, x 2 , x 3 , y, xy, x 2 y, y 2 , y 3 , xy 2 , y 4 , xy 3 , x 2 y 2 , x 3 y, x 4 , y 5 , xy 4 , x 2 y 3 , 
x y ,x y\, 

Q ={l,x, Ax 2 + 3x, 2x 3 + x 2 + Ax, 3y + A, 3xy + Ax + Ay + 3, 

2x 2 y + 5x 2 + xy + 6x + Ay + 3, 6y 2 + 1, 2y 3 + by, xy 2 + 6y 2 + 6x + 1, 
y 4 + 3y 3 + 6y 2 + Ay, 5xy 3 + by 4 + 3y 3 + 2xy + 2y 2 + Ay, 
6x 2 y 2 + xy 2 + x 2 + 6x, . . .}, 
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G ={y 6 + 3y 5 + 2y 4 + 6y 3 + Ay 2 + 5y, 

xy 5 + x 4 y + 6x 3 y 2 + xV + 5xy 4 + 6y 5 + 6x 4 + 2x 3 y + 6z V + 3xy 3 + 
3y 4 + 6x 3 + 6x 2 y + 2xy 2 + Qy 3 + x 2 + 2xy + 6y 2 + x, 
x 2 y A + x 4 y + 3x 2 y 3 + 3xy 4 + 5y 5 + x A + 6x 3 y + 3x 2 y 2 + 2xy 3 + 4y 4 + 
6x 3 + 4y 3 + 6x 2 + 2xy + 3y 2 + a; + 5y, . . .}. 

In the following, several tables show the timings for the computations of 
BM-problems on sets of distinct random points w.r.t. the term order -<i cx or 
-<tdinicx- The algorithms presented in the paper were implemented on Maple 
12 installed on a laptop with 2 Gb RAM and 1.8 GHz CPU. 

Take the field F 2 3, we have 



#2 


200 


300 


400 


500 


BM 


4.968 s 


15.359 s 


34.609 s 


61.172 s 


SPBM 


1.438 s 


3.766 s 


7.141 s 


7.969 s 



For F37, we have 



#S 300 600 900 1200 

BM 16.265 s 121.766 s 420.219 s 1060.203 s 
SPBM 4.172 s 25.125 s 82.000 s 132.719 s 



For F17, we have 





100 


150 


200 


250 


BM 


0.875 


s 2.421 


s 4.953 s 


8.188 s 


GPBM 


0.797 


s 2.125 


s 4.250 s 


5.641 s 


Preprocessin 


g 0.015 


s 0.094 


s 0.172 s 


0.391 s 


#S'/#S 


0.310 


0.393 


0.430 


0.616 


the field F 2 g, we 


have 








#5 


200 


400 


600 


800 


BM 


5.672 s 


38.063 s 


112.156 s 


235.813 s 


GPBM 


5.562 s 


36.906 s 


105.828 s 


135.609 s 


Preprocessing 


0.046 s 


0.313 s 


1.671 s 


8.125 s 




0.125 


0.178 


0.328 


0.711 
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